classdef cDIVASL < handle
    % CDIVASL cDivasl is a class that for diverse-TI
    %   The members include:
    %       TIs:        different Tis in seconds;
    %       dDivasl:    the cASL obj that contains images for different Tis
    %       dM0:        the M0 image
    %       dModel:     the cModel obj that tells which model used
    %                   (buxton's or some others)
    %       dLS,dCMF,dSP:
    %                   the cPerfusion obj that collects the perfusion
    %                   results, including CBF and AAT
    %    The functions include:
    %       cDIVASL(varargin) -  constructions
    %           it can be either constructed using file path to the
    %           original MRI images, or from the asl images and M0
    %           NOTE: M0 is the average of all first repetitions
    %       initialize(obj,varargin)    
    %           initialize the obj (not available now)
    %       setTemplate(obj,varargin)   
    %           set gray or white matter
    %       getImg()
    %           get image (can be replace by directly refer the members)
    %       callLS()                    
    %           compute the CBF and AAT by LS
    %       callCMF()
    %           compute the CBF and AAT by CMF
    % ---------------------------------------------------------------------
    %               VERSION CONTROL 
    % 08/09/2012 *** comments for the class have been added
    %
    % ---------------------------------------------------------------------
    %   
    %
    % COPYRIGHT RESERVED 2012 @ INRIA/IRISA VISAGE, BY LEI YU
    % See also: cASL, cModel, cPerfusion, cMRI
    
    properties
        TIs = [1200:100:2200]/1000;
        dDivasl = [];
        dM0 = cMRI;
        dModel = cModel();
        
        dLS = cPerfusion();
        
        dCMF = cPerfusion();

        dSP = cPerfusion();
    end
    
    methods
        function obj = cDIVASL(varargin)
            if nargin ~=0
                if length(varargin) == 1 % if we have original MRI images
                    obj.dDivasl = cASL(varargin{:});
                    for i = 1:length(obj.dDivasl)
                        obj.dDivasl(i).pTI2 = obj.TIs(i);
                    end
                elseif length(varargin) > 2 % if we just have asl and m0
                    obj.dDivasl = cASL(varargin{1},varargin{2},varargin{3});
                end
                
                M0 = zeros(64,64,14);
                for i = 1:length(obj.dDivasl)
                    M0(:,:,:) = obj.dDivasl(i).dM0.dMRI(:,:,:,1) + M0;
                end
                
                M0 = M0/length(obj.dDivasl);
                obj.dM0 = obj.dM0.initialMRI(M0);
            end
        end
    end
    
    methods
        % Initialize the object (NOT USEFUL NOW)
        function obj = initialize(obj,varargin)
            % Initialize the object
        end
        
        % Set template files
        function obj = setTemplate(obj,varargin)
            nASL = length(obj.dDivasl);
            for i = 1:nASL
                obj.dDivasl(i) = obj.dDivasl(i).setTemplate(varargin);
            end
        end
        % Get Images
        varargout = getImg(obj,command,x,y,slice,rep,maskType);
    end
    
    
    methods
        % This function is to calculate the CBF using Least Square Method
        % maskType = 'Gray' or 'White'
        % dataAmount = 'all' or ':' or '1:20' (should be string)
        function obj = callLS(obj,maskType,dataAmount,selTi,varargin)
            CBF = zeros(64,64,14);
            DT = zeros(64,64,14);
            
            DATA = zeros(64,64,14,11);
            MASK = zeros(64,64,14);
            if strcmp(dataAmount,'all')
                dataAmount = ':';
            end
            for i = 1:11
                DATA(:,:,:,i) = obj.dDivasl(i).getMean('def',:,:,:,dataAmount);
            end

            
            %%%%%%----------------------------------------------------
            %%%%%% Set the M0 and Mask
            M0(:,:,:) = obj.dM0.dMRI(:,:,:,1);
            if isempty(maskType)
                maskType = 'Gray';
            end
            switch lower(maskType)
                case 'white'
                    MASK(:,:,:) = obj.dDivasl(1).dWhite.dMRI(:,:,:,1);% Change the template here
                case 'gray'
                    MASK(:,:,:) = obj.dDivasl(1).dGray.dMRI(:,:,:,1);% Change the template here
            end
            %%%%%%-----------------------------------------------------
            
            if isempty(selTi)
                selTi = 1:11;
            end
            
            
            for s = 3:13
                [h,w] = find(MASK(:,:,s)>0.5);
                selTIs = obj.TIs(selTi)+(s-1)*0.045;
                for i = 1:length(h)
                    dM = DATA(h(i),w(i),s,selTi);
                    obj.dModel.pM0 = M0(h(i),w(i),s);
                    [EST,RESNORM,RESIDUAL,EXITFLAG,OUTPUT] = ...
                        lsqcurvefit(@(x,t) obj.dModel.curvefitting_asl_model(x,t),...
                        [0.0,0.0],...
                        selTIs,...
                        dM(:),...
                        [0,0.001],[300,3.5]);
                    CBF(h(i),w(i),s) = EST(1); %%% Times 3000 (unit mL(100g(-1)min(-1))
                    DT(h(i),w(i),s) = EST(2);
                end
            end
            obj.dCMF = obj.dCMF.initialize(CBF,DT,maskType,dataAmount,selTi,'ls');
            
        end
        
        % This function is to calculate the CBF using Compressive Matched
        % Filter Method
        % maskType = 'Gray' or 'White'
        % dataAmount = 'all' or ':' or '1:20' (should be string)
        function obj = callCMF(obj,maskType,dataAmount,selTi,varargin)
            CBF = zeros(64,64,14);
            DT = zeros(64,64,14);
            
            DATA = zeros(64,64,14,11);
            MASK = zeros(64,64,14);
            if strcmp(dataAmount,'all')
                dataAmount = ':';
            end
            for i = 1:11
                DATA(:,:,:,i) = obj.dDivasl(i).getMean('def',:,:,:,dataAmount);
            end

            %%%%%%----------------------------------------------------
            %%%%%% Set the M0 and Mask
            M0(:,:,:) = obj.dM0.dMRI(:,:,:,1);
            if isempty(maskType)
                maskType = 'Gray';
            end
            switch lower(maskType)
                case 'white'
                    MASK(:,:,:) = obj.dDivasl(1).dWhite.dMRI(:,:,:,1);% Change the template here
                case 'gray'
                    MASK(:,:,:) = obj.dDivasl(1).dGray.dMRI(:,:,:,1);% Change the template here
            end
            %%%%%%-----------------------------------------------------
            
            %%%%%%-----------------------------------------------------
            %%%%% Set the number of TIs
            if isempty(selTi)
                selTi = 1:11;
            end
            %%%%%%-----------------------------------------------------
            
            
            %%%%%%---------------------------------------------------
            %%%%%% Parameter setting for CMF
            deltat = [0.01:0.001:3];
            f = 1;
            %             t = [0:0.01:2.2];
            
            for s = 3:13
                obj.dModel.pSlice = s;
                obj.dModel.pM0 = [];
                [D,AmplitudeD] = obj.dModel.cmf_asl_model(f,deltat,obj.TIs(selTi)+(s-1)*0.045);
                [h,w] = find(MASK(:,:,s)>0.5);
                for i = 1:length(h)
                    dM = DATA(h(i),w(i),s,selTi);
                    obj.dModel.pM0 = M0(h(i),w(i),s);
                    
                    [EST_CBF,EST_DT] = obj.calCBF_CMF(dM(:),...
                        eye(length(dM)),...
                        D*obj.dModel.pM0,...
                        AmplitudeD*obj.dModel.pM0,...
                        deltat);
                    
%                     CBF(h(i),w(i),s) = EST_CBF; %%% Times 3000 (unit mL(100g(-1)min(-1))
%                     DT(h(i),w(i),s) = EST_DT;
                end
            end
            obj.dCMF = obj.dCMF.initialize(CBF,DT,maskType,dataAmount,selTi,'cmf');
        end
        
        
    end
    

    
    
end

